## Figure A10: Effects of Neighborhood Race and Income on Expected Wait Times: City-by- City

## install.packages(c("tidyverse", "fixest"))
## library(tidyverse)

## SET WORKING DIRECTORY HERE
## setwd()

## Loading data
## load("dta.RData")

## Figure A11a (Baltimore)
## Race across-service
race_expected_across <- feols(log_expected_time ~ factor(white_third) |
                                        open_month^open_year, data = dta %>%
                                        filter(city == "Baltimore"), cluster = "geo")
race_expected_across_coefs = as.data.frame(race_expected_across[["coeftable"]]) %>% 
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_expected_within <- feols(log_expected_time ~ factor(white_third) |
                                        service^ open_month^open_year, 
                                      data = dta %>% filter(city == "Baltimore"), cluster = "geo")
race_expected_within_coefs = as.data.frame(race_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Baltimore = race_expected_within_coefs %>%
  bind_rows(., race_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service
inc_expected_across <- feols(log_expected_time ~ factor(inc_third) |
                                       open_month^open_year, data = dta %>%
                                       filter(city == "Baltimore"), cluster = "geo")
inc_expected_across_coefs = as.data.frame(inc_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_expected_within <- feols(log_expected_time ~ factor(inc_third) |
                                       service^ open_month^open_year, 
                                     data = dta %>% filter(city == "Baltimore"), cluster = "geo")
inc_expected_within_coefs = as.data.frame(inc_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Baltimore = inc_expected_within_coefs %>%
  bind_rows(., inc_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A11b (Boston)
## Race across-service
race_expected_across <- feols(log_expected_time ~ factor(white_third) |
                                        open_month^open_year, data = dta %>%
                                        filter(city == "Boston"), cluster = "geo")
race_expected_across_coefs = as.data.frame(race_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_expected_within <- feols(log_expected_time ~ factor(white_third) |
                                        service^ open_month^open_year, 
                                      data = dta %>% filter(city == "Boston"), cluster = "geo")
race_expected_within_coefs = as.data.frame(race_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Boston = race_expected_within_coefs %>%
  bind_rows(., race_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service
inc_expected_across <- feols(log_expected_time ~ factor(inc_third) |
                                       open_month^open_year, data = dta %>%
                                       filter(city == "Boston"), cluster = "geo")
inc_expected_across_coefs = as.data.frame(inc_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_expected_within <- feols(log_expected_time ~ factor(inc_third) |
                                       service^ open_month^open_year, 
                                     data = dta %>% filter(city == "Boston"), cluster = "geo")
inc_expected_within_coefs = as.data.frame(inc_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Boston = inc_expected_within_coefs %>%
  bind_rows(., inc_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A11c (Dallas)
## Race across-service
race_expected_across <- feols(log_expected_time ~ factor(white_third) |
                                        open_month^open_year, data = dta %>%
                                        filter(city == "Dallas"), cluster = "geo")
race_expected_across_coefs = as.data.frame(race_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_expected_within <- feols(log_expected_time ~ factor(white_third) |
                                        service^ open_month^open_year, 
                                      data = dta %>% filter(city == "Dallas"), cluster = "geo")
race_expected_within_coefs = as.data.frame(race_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Dallas = race_expected_within_coefs %>%
  bind_rows(., race_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service
inc_expected_across <- feols(log_expected_time ~ factor(inc_third) |
                                       open_month^open_year, data = dta %>%
                                       filter(city == "Dallas"), cluster = "geo")
inc_expected_across_coefs = as.data.frame(inc_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_expected_within <- feols(log_expected_time ~ factor(inc_third) |
                                       service^ open_month^open_year, 
                                     data = dta %>% filter(city == "Dallas"), cluster = "geo")
inc_expected_within_coefs = as.data.frame(inc_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Dallas = inc_expected_within_coefs %>%
  bind_rows(., inc_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Figure A11d (Houston)
## Race across-service
race_expected_across <- feols(log_expected_time ~ factor(white_third) |
                                        open_month^open_year, data = dta %>%
                                        filter(city == "Houston"), cluster = "geo")
race_expected_across_coefs = as.data.frame(race_expected_across[["coeftable"]]) %>% 
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_expected_within <- feols(log_expected_time ~ factor(white_third) |
                                        service^ open_month^open_year, 
                                      data = dta %>% filter(city == "Houston"), cluster = "geo")
race_expected_within_coefs = as.data.frame(race_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## Plot
race_coefs_Houston = race_expected_within_coefs %>%
  bind_rows(., race_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 



## Income across-service
inc_expected_across <- feols(log_expected_time ~ factor(inc_third) |
                                       open_month^open_year, data = dta %>%
                                       filter(city == "Houston"), cluster = "geo")
inc_expected_across_coefs = as.data.frame(inc_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_expected_within <- feols(log_expected_time ~ factor(inc_third) |
                                       service^ open_month^open_year, 
                                     data = dta %>% filter(city == "Houston"), cluster = "geo")
inc_expected_within_coefs = as.data.frame(inc_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Houston = inc_expected_within_coefs %>%
  bind_rows(., inc_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 


## Figure A11e (Philadelphia)
## Race across-service
race_expected_across <- feols(log_expected_time ~ factor(white_third) |
                                        open_month^open_year, data = dta %>%
                                        filter(city == "Philadelphia"), cluster = "geo")
race_expected_across_coefs = as.data.frame(race_expected_across[["coeftable"]]) %>% 
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_expected_within <- feols(log_expected_time ~ factor(white_third) |
                                        service^ open_month^open_year, 
                                      data = dta %>% filter(city == "Philadelphia"), cluster = "geo")
race_expected_within_coefs = as.data.frame(race_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_Philly = race_expected_within_coefs %>%
  bind_rows(., race_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 


## Income across-service
inc_expected_across <- feols(log_expected_time ~ factor(inc_third) |
                                       open_month^open_year, data = dta %>%
                                       filter(city == "Philadelphia"), cluster = "geo")
inc_expected_across_coefs = as.data.frame(inc_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_expected_within <- feols(log_expected_time ~ factor(inc_third) |
                                       service^ open_month^open_year, 
                                     data = dta %>% filter(city == "Philadelphia"), cluster = "geo")
inc_expected_within_coefs = as.data.frame(inc_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_Philly = inc_expected_within_coefs %>%
  bind_rows(., inc_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 


## Figure A11f (Washington, DC)
## Race across-service
race_expected_across <- feols(log_expected_time ~ factor(white_third) |
                                        open_month^open_year, data = dta %>%
                                        filter(city == "Washington, DC"), cluster = "geo")
race_expected_across_coefs = as.data.frame(race_expected_across[["coeftable"]]) %>% 
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Race within-service
race_expected_within <- feols(log_expected_time ~ factor(white_third) |
                                        service^ open_month^open_year, 
                                      data = dta %>% filter(city == "Washington, DC"), cluster = "geo")
race_expected_within_coefs = as.data.frame(race_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(white_third)3" ~ "3",
                          term == "factor(white_third)2" ~ "2",
                          term == "factor(white_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
race_coefs_DC = race_expected_within_coefs %>%
  bind_rows(., race_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "% White Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 


## Income across-service
inc_expected_across <- feols(log_expected_time ~ factor(inc_third) |
                                       open_month^open_year, data = dta %>%
                                       filter(city == "Washington, DC"), cluster = "geo")
inc_expected_across_coefs = as.data.frame(inc_expected_across[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "across")

## Income within-service
inc_expected_within <- feols(log_expected_time ~ factor(inc_third) |
                                       service^ open_month^open_year, 
                                     data = dta %>% filter(city == "Washington, DC"), cluster = "geo")
inc_expected_within_coefs = as.data.frame(inc_expected_within[["coeftable"]]) %>%
  rownames_to_column(., "term") %>%
  rename(estimate = Estimate, std.error = `Std. Error`) %>%
  select(estimate, std.error, term) %>%
  mutate(term = case_when(term == "factor(inc_third)3" ~ "3",
                          term == "factor(inc_third)2" ~ "2",
                          term == "factor(inc_third)1" ~ "1"),
         term = fct_rev(as.factor(as.character(term))),
         model = "within")

## PLOT
inc_coefs_DC = inc_expected_within_coefs %>%
  bind_rows(., inc_expected_across_coefs) %>%
  ggplot(., aes(x = factor(term), y = estimate*100, group = model, color = model)) +
  geom_pointrange(aes(ymin = (estimate - (std.error * 1.96))*100, 
                      ymax = (estimate + (std.error * 1.96))*100), 
                  position = position_dodge(width = 0.8)) +
  scale_x_discrete(labels = c("Middle", "Bottom")) +
  geom_hline(yintercept = 0, color = "red", linetype = "dashed") +
  scale_color_manual(values=c("darkgray", "black"), name = "",
                     labels = c("Across Service", "Within Service")) +
  labs(x = "Per Capita Income Tercile",
       y = "% Change in Expected \n Wait Time vs. Top Tercile") +
  theme_classic() +
  theme(panel.grid.major = element_blank(),          panel.grid.minor = element_blank(),         panel.background = element_rect(colour = "black"),         legend.position = "bottom") 
